ON CONVERGENCE OF THE PROJECTIVE INTEGRATION 
METHOD FOR STIFF ORDINARY DIFFERENTIAL EQUATIONS 
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Abstract. Wc present a convergence proof of the projective integration method for a class of 
deterministic multi-dimensional multi-scale systems which are amenable to centre manifold theory. 
- - - The error is shown to contain contributions associated with the numerical accuracy of the microsolver, 

m : the numerical accuracy of the macrosolver and the distance from the centre manifold caused by the 

1 combined effect of micro- and macrosolvers, respectively. We corroborate our results by numerical 

' simulations. 
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' 1. Introduction 

i Devising efRcient computational methods to simulate high-dimensional complex 

systems is of paramount importance to a wide range of scientific fields, ranging from 
nanotechnology, biomolecular dynamics and material science to climate science. The 
dynamics of large complex systems is made complicated by their high dimensionality 
and by the possible existence of active entangled processes running on temporal scales 
' varying by several orders of magnitude. To resolve all variables in a high-dimensional 

. system and capture the whole range of temporal scales is impossible given current 

computing power. In many applications, however, the modeller is only interested in 
the dynamics of a few relevant slow macroscopic coarse-grained variables. How to 
extract from a dynamical system the relevant dynamics of the slow degrees of free- 
dom while ensuring that the collective effect of the unresolved variables is implicitly 
represented is one of the most challenging problems in computational modelling. 
There are two separate scenarios when such a dimension reduction is possible: 
scale separation and weak coupling We restrict ourselves here to the big class 



in 
oo 



' of time scale separated systems, in particular to deterministic stiff dissipative systems. 

. Recently two numerical methods to deal with multi-scale systems have received 

CO ' much attention, the projective integration method (PI) within the equation- free 

framework [1, [3l and the heterogeneous multiscale method (HMM) 0, 0]. These 
powerful methods have successfully been applied to a wide range of problems, 
. including modelling of water in nanotubes, micelle formation, chemical kinetics and 

'jjj ' climate modelling [isl . H, [l^l . The general strategy of both methods is the following: 

. Provided there exist closed (but possibly unknown) equations for the slow variables, 

the simulation is split between a macrosolver employing a large integration time-step 
and a microsolver, in which the full high dimensional system is integrated for a short 
burst employing a small integration time step. The two methods differ in the way 
the information of the microsolver is utilized in the macrosolver to evolve the slow 
variables. In the equation-free approach this is done without any assumptions on the 
actual form of the reduced equations, using the microsolver to estimate the temporal 
derivative of the slow variable, which is then subsequently propagated over a large 
time step. The underlying assumption is that the fast variables quickly relax to 
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their equilibrium value (conditioned on the slow variables), and that subsequently 
the dynamics evolves on a reduced slow manifold. In heterogeneous multiscale 
methods, on the other hand, information available (i.e. through perturbation theory 
such as for example centre manifold theory, averaging theorems or homogcnisation 
[II, m , 15 , 1^ ) is used to determine the functional form of the reduced slow vector field; 



the microsolver is then used to estimate the coefficients of this slow vector field con- 
ditioned on the slow variables. For further details we refer the reader to [l2l. [isl lil [iij . 

The question we address here is the convergence of the numerical PI approximation 
to the true solution of a full multi-scale system of stiff deterministic ordinary 
differential equations. There exists a body of analytical work on the convergence of 
HMM 0, 0, including formulations for stochastic differential equations 0, flTj . 
Convergence results for PI were obtained for specific deterministic systems and for 
on-the-fly local error estimates in numerical algorithms [1, [l^ . Stochastic differential 
equations were treated in [l^l- Here we provide global error estimates for PI for a 
subclass of systems amenable to centre-manifold theory. 

The paper is organized as follows. In Section [2] we discuss the class of dynamical 
systems studied, and briefly summarize in Section [3] HMM and PI for these systems. 
In Section m the main part of this work, we derive rigorous error bounds for those nu- 
merical multiscale methods. In Section [5] these are numerically verified. We conclude 
with a summary in Section [51 

2. Model 

We consider deterministic multiscale systems for slow variables £ M" and fast 
variables S M™ of the form 

ye = g{xe,ye), (2.1) 

i, = i(-Aa-, + /(y,)), (2.2) 

with time scale separation parameter < e <C 1. Without loss of generality we assume 
that {xcUe) ~ (0,0) is a fixed point. We assume there is a coordinate system such that 
the matrix AsM™^™ is diagonal with diagonal entries Xii>0. We further assume 
that we allowed for a scaling of time such that min(Aii) = 1 and define max(Ai,;) = A. 
We assume that the vectorfield of the slow variable g{x,y) is purely nonlinear with 
detDg(0,0) ^ 0, where Dg denotes the Jacobian of g{x,y), so that there exists a centre 
manifold x = h^{y) on which the slow dynamics evolves according to 

r = G(r), (2.3) 

with Y ^y + 0{e). The reduced slow vectorfield is given by (see for example [l[) 

Giy)^gihM,y)- (2.4) 

Initial conditions close to the fixed point are exponentially quickly attracted towards 
the centre manifold along the stable manifold of the fast variables x^- Near the 
centre manifold the dynamics slows down and is approximately determined by the 
dynamics of the slow variables ()2.3p only. Centre manifold theory assures that on 
times T^O{l) the dynamics of ()2.ip is well approximated by the reduced slow 
system (|2.3p (see for example [H). 
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It is "worthwhile to formulate centre manifold theory in the framework of averaging 
, and view the effect of the fast variables on the slow variables through their induced 
empirical measure which is approximated by ^y{dx) = 5{x — h^{y))dx , conditioned on 
the slow variables. The dynamical system (|2.ip - (j2.2p can equivalently be described by 
its associated Liouville equation for the probability density function p{x^,y^). After 
an initial transient, the solution of the Liouville equation is approximated by 

p{x, y) = 5{x - h, {y))p{y) + 0{e) . (2.5) 

We can now reformulate (|2.4p as 



G{y) = j g{x,y)py{dx) 

^g{K{y),y). (2.6) 

The underlying assumption is that the measure p,y{dx) = 6(x — h^{y))dx is the physi- 
cally observable measure on the y-fibre; that is, Lebesgue almost all initial conditions 
of the fast variables will evolve to generate fiy{dx), if sufficiently close to the centre 
manifold. Centre manifold theory [l[ makes this approach rigorous and p{y) is the 
invariant density of the reduced system (|2.3p . This is exploited explicitly in HMM 
and implicitly in PI. In the following we will approximate the centre manifold h^[y) 
hyK{y)^K-^f{y) + 0{e). 

3. Numerical Multiscale Methods 

We consider here two numerical methods to deal with the deterministic multi- 
scale system (|2.ip - (|2.2p . namely heterogeneous multiscale methods and projective 
integration. We use the formulation of PI as in [s, 13| and of heterogeneous multiscale 



methods as in [1, 0] ■ We will see that both methods can be formulated in the same 
framework. 



We denote by x" and j/" the numerical approximations to the solutions of (|2.ip - 
(|2.2[) evaluated at the discrete time f", ^^(i") and ydt"') respectively. Further 
we denote by Y(t") the time-continuous solution of the reduced ordinary differ- 
ential equation p.3|) evaluated at time t". Throughout the paper superscripts 
denote discrete variables whereas brackets are reserved for continuous variables. 
We will employ a forward Euler method for both the micro- as well as the macrosolver. 



3.1. Heterogeneous Multiscale Methods 

The microsolver consists of M microsteps with micro-time step St 0]. The slow 
variables are held fixed during the microsteps, assuming infinite time scale separation. 
Let us denote by a;"'™ the mth microstep of the fast variables which was initialized at 
the nth macrostep at = nAt with = j;""^'^-'^ = and = y". The microstep 
for a forward Euler scheme then is 

^n,m+l ^ ^n,n, _ ^ (^.t"'™ - /(y")) , 

for m = Q,...,M. Invoking the Birkhoff ergodic theorem, the time series of the fast 
variables x""™ is used to estimate the average in the reduced slow veetorfield (|2.6p . 
and the macrosolver is initialized at the previous macrostep y" and is written for a 
forward Euler step as 



y"+i=j/" + Ai5(x",y") 



(3.1) 
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with the time-discrete approxiraation of the reduced slow vectorfield 

M 
m=0 

for some weights Wm- This is justified as long as the slow variable remains constant 
while integrating over the fast fibres. Being associated with the empirical approxima- 
tion of the invariant density induced by the fast dynamics (cf. (|2.6|) '). the weights Wm 
satisfy the normalization constraint 

M 

Y,Wra = l. (3.2) 

m=0 

The weights should be chosen to best approximate the invariant density of the fast 
variable. A natural choice for our system (|2.ip - (|2.2p . where the fast dynamics rapidly 
settles on the centre manifold, and does not advance significantly on the centre man- 
ifold in M microsteps is to set 

{0 for TO < M 
1 tor TO = M 

The convergence of this scheme has been analyzed in Q . More general weights were 
discussed in Q. 

It is pertinent to mention that for this formulation of HMM the fast variables have 
to be initialized before each application of the microsolver; see for example Q for 
a discussion on initializations. In PI, as described in the next subsection, the fast 
variables are initialized on the fiy. 

The seamless heterogeneous multiscale method 0, 0] can be applied to the case when 
the slow and fast variables are not explicitly known. The seamless HMM formulation 
for our system ()2.1|) - (j2.2|) advances both, fast and slow, variables simultaneously first 
through M microsteps, and then subsequently through one macrostep. We adopt 
here the formulation described in Q for deterministic systems. Introducing y"^™ 
analogously to cc"'™ the M microsteps are executed as 

^.n.m+l ^ ^n,m _ ^ (Ax"''" - /(?/"'™)) (3.4) 

y"'™+i = y"'" + (S<5(a;"'™,y"'™). (3.5) 
The macrostep is initialized at (a;",y") = (a;"'°,i/"'°) and takes the form 



X 



At 

= V Wm (Ax"'™ - /(y"''")) , (3.6) 

e ^-^ 

m=0 

y"+i = y" + Aig(x",y"), (3.7) 
with the time-discrete approximation of the slow vectorfield 

g(x",y") := M^„.5(x"'™,y"'™) . (3.8) 

m=0 

Again, the weights Wm need to satisfy (|3.2p (see for choices of weights Wm)- 
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3.2. Projective Integration Wc present a formulation of projective integra- 
tion which is close to the seamless formulation of HMM. PI advances both slow and 
fast variables at the same time, as done in seamless HMM. However, the PI scheme 
proposed in 0] does not start the macrostep at (a;"'°,2/"'°) as done in our adopted 
formulation of seamless HMM, but at (a;"'*'',?/"'*^), allowing for nontrivial dynamics 
of the slow variables during the M microsteps. This takes into account finite time 
scale separation e, ignored in HMM. The microsteps of PI are exactly as in seamless 
HMM, (|3.4p - (|3.5p . The macrostep of PI is given by estimating the temporal time 
derivative of the slow (and fast) variables using the microsolver according to 



j.n.M _ n,M-l 

-At ; 

6t 



„.n,M _„,n.M-l 

y»+i = y'M/ + Ai^ 1 , 

which becomes upon using the Euler step of the microsolver p.4p - p.5p 

^n+l ^ ^n,M „ ^(^^.rM/ „ /(y".^^)) (3.9) 
yn+l^ynM_^^^g(^^n,My^,M-^^ (3.10) 

Here the time between two macrosteps is ^At + AIdt and wc have t" ^ntj^. 
Upon substituting the microsteps p.4p - p.5p into the macrosolver (|3.9p - (|3.10p we 
reformulate PI, such that it formally resembles seamless HMM, as 

M 
^ m=0 

y"+i=y" + t^5(x",y"), (3.12) 
with the weights Wm now defined as 

St r 

— tor m < M 

^^-^{^It ■ (3-13) 

— for ni = M 
tA 

Hence PI with a microstep of dt and macrostep of At is equivalent to seamless HMM 
with a microstep of St, a macrostep of tA and a particular choice of weights Wm- 

4. Error analysis for Projective Integration 

We will now provide rigorous error bounds for PI in the formulation p.lip - p.l3p . 
We follow the general line of proof used in for error bounds of HMM with the 
choice of weights ()3.3p . 

Throughout this work we assume the following conditions on the global growth and 
smoothness of solutions of our system and on the numerical discretization parameters 
of PL 

Assumptions 

Al: The vectorfield /(y) is Lipschitz continuous; that is there exists a constant Lf 
such that 

\.f{yi)- f{y2)\<Lf\yi-y2\ ■ 
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A2: The vectorficld g{x,y) is Lipschitz continuous; that is there exists a con- 



stant Lg such that 

\9{xi,yi) ~ g{x2,y2) \ < Lg{\xi - x2\ + \yi - y-il) ■ 

A3: The vectorfield f{y) is bounded for all y; that is there exists a constant Cf 
such that 

Cf=sup\f{y) \ . 

A^: The vectorfield g{x,y) is bounded for all x,y] that is there exists a constant Cg 
such that 

Cg = sup|g(x,y)| . 

A5: The reduced vectorfield G{Y) of the slow dynamics is of class C^; that is 
there exists a constant C* such that 

C*=sup|f(t)| . 

A6: The microstep size 5t resolves the fastest of the fast variables, while the 
niacrostep size At resolves the dynamics of the slow system, so that 

2e 

< (5t < — < At . 
A 

A 7: The total time M5t of the microsteps is sufficiently short so that 

LcMSt < Lg{l + Lf)M5t < 1 , 
where Lq < Lg{l + Lf) is the Lipschitz constant of the reduced dynamics 



A8: The macrostep size At, number of microsteps M and microstep size St are 
chosen such that 

/ MSt\ e .r ^ c 2e 
Atexp <- if 0<St<' 



e 



A -A + 1' 



or 



f M6t*\ e 26 ^ 2e 

At exp < - if - — - <St< — . 

V e J X A+1 A 

with 6t* = 2e-XSt. The range 2e/(A + 1) < (5t < 2e/A corresponds to 
< 5t* <2e/{X + l). This assumption is necessary to bound the distance of 
the fast variables from the centre manifold over the macrosteps. 

The global Lipschitz conditions can be relaxed to local Lipschitz conditions by the 
usual means. 

We will establish bounds for the error between the PI estimate y" and the solution 
of the full system 2/e(t") 
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Our main result is provided by the following Theorem: 

Theorem 4.1. Given the assumptions (Al)-(A8), there exists a constant C such 
that on a fixed time interval T , for each n such that nt/^ < T, the error between the 
PI estimate and the exact solution of the full multiscale system \2.1\) - ^K^] is given by 



tA / 7 "A + 1 



where d" measures the maximal distance of the fast variables x from the approximate 
centre manifold x^ f{y) := A^^/(y) over all macrosteps and is estimated by 



\d"\< 



with 6t* ^2e — XSt. 



s^AtXe-^ -A + 1 

e-AtXe — — A + 1 A 



Before proceeding to the proof of the theorem, it is worthwhile to interpret the 
bounds on f The term proportional to t^ = At + MSt reflects the first order conver- 
gence of the forward Euler numerical scheme used to propagate the micro- and macro- 
solver respectively. The terms proportional to the time scale parameter e represent 
the error made by the reduction as well as an additional error incurred during the drift 
of the slow variable over the microsteps. The term proportional to exp(— A/(5t/£)|(i"| 
measures the exponential decay of the fast variables towards the slow manifold. 

4.1. Error Analysis We split the error into two parts according to 

5" = |ye(t")-y"| 

<\y,{t")-Yit")\ + \y"-Yit-)\, 

where the first term describes the error between the exact solution of the full system 
(|2.ip - p.2p and the reduced slow system p.3p . which we label reduction error, and 
the second term the error between PI and the exact solution of the reduced slow 
system ()2.3p . which we denote by discretization error. We will bound the two terms 
separately in the following. 

4.2. Reduction error Defining the error between the exact solution of the 
full system p.ip - p.2p and the reduced slow system p.3p 



\E-\^\y,{f^)-Yr)\, (4.1) 

and setting the initial conditions close to the slow manifold with i/e(0) ^Y{0) + co^y£ 
and x^(0) ~ f {y^(0)) + cq^x , we can formulate the following theorem: 

Theorem 4.2. Given the assumptions (A1)-(A3), there exists a constant C'l such 
that on a fixed time interval T , for each < T , the error between exact solutions of 
the reduced and the full system is bounded by 



\E:\<Cie, 
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with 

Ci=niax(|2/,(0)-r(0)|,L/L<,C/",i<,|de(0)|))e^«(i+^^)*", 

where de = Xg ~ fiVe) measures the distance of the fast variables from the approximate 
slow manifold. 

Proof. The proof is standard in centre manifold theory and is included for com- 
pleteness. Using assumptions {A1)~{A2) on the Lipschitz continuity of / and g we 
estimate 

j^\y,{t)-Y{t)\ = \g{x,,y,)^g{f{Y),Y)\ 

<\g{xe.ye)~9{^e.Y)\ + \g{x,,Y)-g{f{Y),Y)\ 

<L,\y,^Y\+Lg\x,-f{Y)\ 

<Lg\y,-Y\+Lg\x,- f{y,)\ + Lg\f{y,) ~ f{Y) \ 

<Lg{l + Lf) \y, -Y\+Lg\x, - f{y,)\ , (4.2) 

To estimate the second term, we differentiate the distance of the fast variable to the 
slow manifold = — f{ye) with respect to time, using (|2.2p . as 

de ^ -jde - Df{y^) g{xe,ye) , 

where D/ denotes the Jacobian matrix of /. Integrating and using the boundedness 
assumptions (A^) on g and the Lipschitz continuity assumption (Al) on /, we readily 
estimate 



|4(OI<exp(--t)Me(0)| + 



A 



exp( {t-s))Df{y=:{s))g{xeis),yeis))ds 

£ 

<e"^ 14(0)1+ sup \Dfg\ I e~'-^ ds 

\g\<Cg Jo 

<e-i\deiO)\+sLfCg, (4.3) 

which signifies the exponential attraction of the fast dynamics towards the slow centre 
manifold. Substituting (|4.3p into (|4.2p yields upon application of the Gronwall lemma 



\yeit)~Yit)\ < e^«(i+^/)* (|ye(0) -F(0)| +eL^L3C3i + £Lg|4(0)|) , 

which immediately yields the desired result using the initial condition de{0) = cq_x and 
y,(0) = r(0)+co,,£. □ 

4.3. Discretization error In order to bound the discretization error |?/" — 
F(i")| we first estimate how close the fast variables remain to the centre manifold 
during the application of the microsolver. Defining the deviation of the PI approxi- 
mation of the fast variables from the approximate centre manifold at time i" + mdt 

d"'" = x"'"-/(y"'™), (4.4) 

we formulate the following Lemma which constitutes the discrete time version of ( 
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f{y) 



X 



,n,m— 1 „ n,m — 1 



) 



y 



Fig. 4.1: Illustration of equation (|4.5p for XeGM, so that A = l. The two dots show 
one microstep in PI. 



Lemma 4.3. Given assumptions (Al), (A4-) and (A6), the error between the fast vari- 
ables and the approximate centre manifold during the application of the PI microsolver 
is bounded for all <m< M by 



The first term is a manifestation of the exponential decay of the fast variables towards 
the slow centre manifold along their stable eigendirection. The second term propor- 
tional to e describes, as we will see below, the cumulative effect of changes of the slow 
variables y during the microsteps causing a departure from the centre manifold for 
nonconstant /(y). 

Proof. Using the definition of the PI microsolver (|3.4p . a recursive relationship 
for d"'™ is readily established as 



Let us pause here for a moment to discuss the terms appearing in this equation. The 
first term illustrates that at each microstep the distance between the ith fast variable 
and the approximate slow manifold decreases by a factor of l — StXa/s. The second 
term represents the change in the approximate slow manifold during the drift of the 
slow variable y in one microstep. This is illustrated in Figure HTT] 




— X' 




(4.5) 
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The cumulative effect of A'l mierosteps on d"-™ is easily obtained from (|4.5p as 

m—l / f, \ k 



St 



)-Kf' 



)) 



(4.6) 



Taking absolute values and using the Lipschitz condition on / (Al), the bound on 
the slow vectorfield g (A4) and p.Sp . the cumulative effect of M mierosteps on the 
distance of the fast variables to the approximate centre manifold is 



< 



e 

£ 



-LfCgSt 



-LfCgSt- 



m — l 

E 

k=a 
1-1 



St . 
A 

£ 

St \^ I m 



1_|I_^A| 

To avoid divergence from the centre manifold we require 



. St . 
I A 

e 



, St ^ St , 
- max< 1 , A 1 > < 1 . 

£ e 



implying the standard condition < (5t < 2e/A. For < (5t < 
we obtain the bound 



< 1 



St_ 

e 



Tl,0| 



-eLfCg 



2e 
A+1 ■ 



St 

6 



-#aUi- 



st 



and 



(4.7) 



Furthermore for -^<St<^, introducing (S<* = 2£ - A(5t, |l- ^A| = A^ - 1 = 1 - 

A+l A ' ° ' I e I £ e ' 

yielding 



< 1 



St* 
e 



\d 



St 



St 



1 



1 



St* 

e 



(4.8) 



Lemma follows from (|4.7p - (|4.8p . Note that the term proportional to e was gener- 
ated by a sum of error terms in St relating to the change in the manifold during the 
drift of the slow variable y"'™ over each microstep. □ 



Remark 4.1. In the limit St* -^0 (i.e. St^2e/X) the term 

I, which is neglected in Lemma \4.3\ is crucial to obtain a finite estimate 
mStL f Cg . 



ni- 

in 

"n,m\ ^ 



Remark 4.2. In the case gM with A= 1, a direct estimation of |^.5p implies that 
for the particular choice St = e we have < LfCgSt, i.e. there is no dependence on 

initial conditions for the fast variable; this is a peculiarity of the underlying forward 
Euler scheme used (cf. \3.4^ ) where for St~e we have — /(y"'™). 

We have bounded the distance of the PI approximation of the fast variables x"'™ 
from the slow centre manifold /(j/"'™) over the mierosteps. We now establish bounds 
on the distance of the PI approximation of the slow variables y"'"* from the solution 
of the reduced dynamics on the centre manifold over the mierosteps. We therefore 
introduce the auxiliary time-continuous system 



y(")(s)-G(y("'(s)) 



(4.9) 
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(4.10) 



which resolves the reduced slow dynamics (j2.3p initialised after each macrostep at time 
t" ~nt^ with the PI approximation of the slow variables y" . Further we introduce 
the discrete Euler approximation of the reduced slow dynamics (|2.3p 



(4.11) 
(4.12) 



for arbitrary initial conditions rj. Note that for the particular choice rj = y" (j4.1ip - 
(|4.12p corresponds to an Euler discretization of the system (|4.9p - (|4.10p . 



Lemma 4.4. Assuming (Al),(A2),(A5) and (A7), the numerical estimate y"^™ of the 
slow variable is close to ip'^n with 



3Lge\d 



in.Ol 



< 



iri.Ol 



■^ + 0{mSt^) 



if 0<(5t< 



2e 
A+T 



■f 2e 2e 
if - — - <St< — 
A + 1 A 



Lf + 1 J St* 

Proof. Employing the auxiliary system (|4.9p - (|4.10p . we split the error into 

- ¥>™ I < - r(") (mSt) \ + |y (mSt) ~ ^™ | . (4.13) 

The second term above can be readily bounded by the standard proof for Euler 
convergence (see for example [HI). The first term is bounded by a slight variation of 
the same proof. For completeness we provide the proof of the bounds for both terms, 
beginning with the bound on the second term. By Taylor expanding we write upon 
using the auxiliary dynamics (I4.9p 



r (mSt) = y {{m-l)St)+ StG{Y^''^ {{m - l)St) ) +0{St^ 



(4.14) 



Note that assumptions (Al) and {A2) imply a Lipschitz constant LG<Lg{l + Lf) 
for the reduced vectorfield G. Using the Euler discretization (|4.1ip of the reduced 
auxiliary system (|4.9p - (|4.10p with initial condition y", and employing assumption 
{A5) with C* =sup|f(t)| to bound the 0((5t^)-term in gH]), we obtain 



I y (") (mji ) - v"" I < y ^"H ("^ - 1 ) '^^) - '/'IT"" ' 

+ 5i|G(y(")((m-l),5f))-G((^;ri) 

< (1 + LaSt) I {im-i)dt)- (^;'r 1 

<C*5t^Y.{l + LGSt)\ 

k=0 

since we initialize with F ^"^ (0) = <^y„ . Evaluating the sum we find 



-C*St^ 



+c*st^- 



|y(")(m(5i)--'y3" 



<2C*m(5r , 



(4.15) 
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where assumption {AT), that LQMSt< 1, was used to bound e^''™** — 1 < 2LQm5t. 
To bound the first term in (|4.13p . |?/"'™ — F'^"'(m(5i)|, we analogously obtain 

+ (S4g(a;"'"'-\y"'"-i)-G(y"'"-i)| 
-5t G(zj"''"-i)-G(r(")((m-l)(5t)) +0*5^ 



< 



(l + LaSt) 2/"'"-i-y(")((m-l)(5i) 



LgStld 



(4.16) 



Employing Lemma [4.31 with 0<St< we obtain 



7n — 1 



i-l-fe 



k=a ^ ^ ' 

rn — 1 

+ [LgLfCgeSt + C*St^] ^{l + LGSt)\ 

k=0 

since = = y(") (0). Evaluating the sums and using assumption (AT) yields 



5t 



+ 2LgLfCgemSt 

< Lge\d''-" 1(1 + 2LG7nSt) + 2C*m6t^ + 2-^^*^^^ 



L/ + 1 



<3L„e|d"'"|+2C*m5r + 2 



2 , o^/'^s^ 



Lf + 1 



which concludes the proof of the lemma for < i5t < . Analogously for -^p^ <6t< 
Lemma is used to estimate (|4.16p 



y(")(m^t) <5tig|d"'°| ^(I + Lg'^O'' ( 1 

fc= 

(5i 



fc=0 
j-2 



6t 

e 



LgLfCge— + C*St 



St* 



7n — 1 



J2{l + LGSt)>- 



-2G*TO(5r 



(5t 

+ 2LgLfCgem— 



□ 

Lemma 14.41 established bounds on the PI approximation of the slow variables 
and an Euler approximation of the reduced system during the microsteps. We will 
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now establish bounds on the PI approximation of the slow variables and an Euler 
approximation of the reduced system during one maerostep. In order to do so we 
introduce the vector field 



M 



G(r7):=^iy„G(^™), 



(4.17) 



m=0 



with weights Wm given by p.l3p as in PI. We now proceed to bound — 
G(y")|. 

Lemma 4.5. Assuming (A1)-(A6), the auxiliary vector field G is close to the vector 
field g of PI, with 



tA 



+ 3L,il + Lf)e)\d--"\ 



-iLgLfCgE 



|g(x",2/")-G(y")|<<^ 



-LgC- 



-\d 



if < St < 



2e 
A+T 



if T <St< — 

A + 1 A 



Proof. Employing the Lipschitz condition (Al) on / and {A2) on g^ we write 



\~g{x\y-)-G{y")\ 



M 



771 — 

M 

m=0 
M 

< E W,^Lg 



Y,W„7 g(x"^™,y"^")-G(^^^) 



ni=0 



M 



<J2w,nLg |d"'"| + (l + L;)|y 



(4.18) 



Using Lemmas 14.31 and 14.41 for the case when Q<5t< , we expand 



M 



\~g{x-,y-)-G{y-)\<J2^" 



6t 



+ 3Ll{l + Lf)e\d 



-iLgLfGgE 



iri,0 



(4.19) 



Inserting the weights W,n p.l3p which characterize PI, and evaluating the geometric 
series, we obtain 



M-l 



5t 



)[x-,y-)-G{y-)\<Y,-Lg[l- 

m— 



5t_ 

e 



At 



Lg[\ 



M 



_7Tl.0| 
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+ 3LgLfCge + 3Ll{l + Lf)e\d'''°\ 
<Lgi— + 3Lg{l + Lf)s + e- — 



|d"^"|+3LgL/Cg£, 

2e ^ ^ 2e 



which concludes the proof for that case. Using Lemmas 14.31 and 14.41 for <St<^ 
in (|4.18p . we obtain analogously 



M 



|5(x",y")-G(j/")|<^Ty„ 



-3^l2(1 + L/)£M"^°| 



'^JjrLgLfCge 



(4.20) 



(5i 



Remark 4.3. The above estimates (|4.19p and ()4.20p not onZj/ /loM /or PI, hut also 
for the seamless formulation of HMM, since they are independent of the particular 
choice for the weights Wm and of the method of reinitialization of the fast variables 
x". 

We have established in Lemma 14.51 that the vector field G given by (|4.17p is close to 
the PI approximation g of the slow dynamics over a time step of tA- In the following 
Lemma wc demonstrate that G can be used to step forward the reduced slow variables 
in a macrostep. 

Lemma 4.6. G(y(t")) provides a numerical estimate of the reduced slow vector field 
with 

y (r+i) ~~ Y{f') = tAGiY(f')) + o{ti) , 

where the error term 0{t\) is bounded by 2C*t\. 
Proof. By Taylor expanding we write 

r(t"+i)-y(r) = tAG(y(r))+<A(G(r(r))-G(y(r)))+o(ti), 

where according to assumption {AS) the C'(iA)-tcrm is bounded by C*t\. We now 
bound 



G(y(r))-G(y(r)) 



M 



M 

<Y,WrnLG\Y{tn-^YiV^) 

< LaGatA , 



where Cq < Cg is the maximum of |G| . Noting that C* ~ sup jF | = sup |_DG(F) G{Y) \ = 
LgGg completes the proof. □ 

Using Lemmas 14.31 14.41 14.51 and 14.61 we now proceed to bound the discretization error 



\E2\ = W'~Y{t^)\ 



(4.21) 
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and formulate the following theorem. 

Theorem 4.7. Given the assumptions (A1)-(A8), there exists a constant C2 such 
that on a fixed time interval T , for each nAt < T , the error between the solution of the 
projective integration scheme and the exact solutions of the reduced system is hounded 
by 



tA J J - \ + l 



where :=maxo<i<„_i|d' 

Remark 4.4. The error estimate involves the well known exponential decay of the 
fast variable towards the approximate centre manifold, leading to a loss of memory of 
the fast initial condition; PI, however, involves an additional error term proportional 
to the maximal distance of the fast variable to the approximate centre manifold which 
involves no memory loss with M -^00. 

Proof. We estimate the difference i?^' = y" — applying Lemma [4.61 

= (G(y"-i) - G(y(r-i))) <A + - G(y"-i)) t^ + 0{tl) 

^Clr^Ell-hA + a.^-itA- (4.22) 
We used the mean value theorem for vector-valued functions to introduce 

C^:^ f DG(r(r) + 6i(?/"-y(t")))d6', 
Jo 

where DG is the Jacobian matrix of G, and we set 

a„:=5(x",y")-G(y") + 0(tA), 
where according to Lemma 14.61 we have 

|a„|<l5(a^",y")-G(y")| + 2G*iA. (4.23) 
Taking absolute values we obtain 

ii?,"i<(i+iAii/:r'ii)i^r'l+i"-ii^A. 

Employing assumption {A4) on the boundcdness of g, we define 

a= max laA , 

0<j<ri-l 

and using {A5) on the boundcdness of G it is easy to show that 

= , max \\Cg\\- 
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We obtain, evaluating the geometric series, 



n-1 



(4.24) 



where we used that at ri = we initiahzc with £'2 = 0. Recahing the bound (|4.23p 
for \a„\ and employing Lemma [4.51 for {)<5t<^^ we obtain for the bound of the 
discretization error 



\E2\<- 



2C*tA+ Lq 



tA 



^3Lg{l + Lf)e + e-— \d"\+3LgLfCge 



(4.25) 



Employing Lemma 14.51 for <6t<^ we analogously obtain 

|i.."|<^{2C*.. + £,(^^ + 3^X,(l + X,)s-,e--|- 



\d" 



(4.26) 



Theorem HJ] follows realizing 6t* < St for St > 2e/(A + 1). □ 

Besides the parameters used in the numerical scheme, i.e. the macrostep size Ai, 
the number of microsteps M with microstep size St, and the time scale parameter e, 
the error bound also involves the maximal deviation of the fast variables to the 
approximate centre manifold. 

The bound on is established in the following Lemma. 

Lemma 4.8. Assuming (Al), (A4), (A6) and (A8), the maximal deviation of 
the fast variables from the approximate centre manifold over n macrosteps \d"'\ = 
maxo<i<„_i Id"'"! satisfies 



if < St < 



|d"|< 



,0,01 , eLfCg{l + X)At 

e-AtXe — — 
,0,01 , eLfCg{l + X§r)At 2s 



2e 
A+T 



AtAe- 



A + 1 A 



Proof. We bound for 0<i<n 

W''\ = W-nyl\ 



i-l,M 



At 



(A: 



<At-\d 

£ ' 



+ \.f{y'-'^^'')~f{yl\ 



fif~''''))-fiyl 



where we have used assumption ( A 6) to simplify 



At^-1 



= At--l<At-. 

£ £ 
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Employing Lemma [4.31 for 0<5t< ^p^, and assumptions {Al) and {A4), we obtain 

•^1 n~l.M\ , r . L,i-1,M „,i| (4 27) 



A M5t 

<At-e 

e 



|d*-i^°|+L/Cg(l + A)At 
Evaluating this recursive relationship yields 

1 



|d*'0|< 



A MSt 

At-e ~ 

e 



LfCg{l + X)At- 



At^e — - 



1-At^e- 



(4.28) 



LfCgil + X)At 
1-At^e-^ 



where we used assumption {A8) with <dt< -^p^- Using Lemma 
in (|4.27p we obtain 



for 



2e 
A+1 



<St<^ 



< 



. A Af at* 

At-e — ~ 



St 



1 



At^e 



< d' 



0,01 



M-|+%C,(l + A^)At-^_^^^^_ 

e 

St 



(4.29) 



e%Cg(l + A#)At 



<5t* 



£ - AtAe 



where we used assumption {A8) with < (5t < □ 

Remark 4.5. In the limits of Assumption (A8), Atexp(— 



At exp 



MSt* 



— >■ y , the terms 1 — 



At^e — — 



At^ 



m ()4.29p , which are neglected in Lemma \4-8[ are crucial to obtain finite estimates 
\d"\< \d°'° I + nAtLfCg (1 + A) and |d" | < \d"-° \ + nAtLfCg (l + 1^ A) , respectively. 

Theorem l4.1l follows by combining Theorems l4.2l and l4.7l with Lemma H?51 realizing 
St*<St for dt>2e/{\ + l). 

5. Numerical confirmation of the error bound \EJ}\ 



We now illustrate the error bound \E2\ of Theorem 14. 7[ (|4.25p . which we recall 
here including the constants as obtained in the proof 



\E2\<3-^lc*tA + Lg^- + Lg{l + Lj)s + e-^ 



\d^^\+LgLfCge\, (5.1) 



for 0<6t<^pj, and 

\E2\<3^^^lc*t^ + Lg^- + Lg{l + Lf)e + e-^ ]\d-\+LgLfCgs[ 



(5.2) 
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for^<<5t<f. 

We demonstrate the scaling of \E2\ with respect to the macrostep size At, the 
timescale parameter s and the reinitialization error of the fast variable for fixed 
final time T = nt^ ■ Note that the dependencies of e and At are complicated through 
the dependency of on those parameters. We show results for simulations using 
the following multiscalc system 

ye = -Xeye-ayl (5.3) 
.^^-x^+^^^by^ ^ (5.4) 

which has, at lowest order in e, the slow limit system 

Y^-Ysw?{hY)~aY'^ . (5.5) 

For higher order approximations of the centre manifold and the associated coordinate 
transformations relating y and Y the reader is referred to the very useful webpage 
[m (see also 0). 

The system l|5.3p - (|5.4p with initial condition y^{Q)>Q is locally Lipschitz with Lips- 
chitz constant Lf<b and Lg = max(|j;e| + 2a|ye|) where the maximum is taken over 
the local region around the fixed point at (xg^ye) = (0,0) under consideration. The 
vector field of the slow dynamics ()5.3p is locally bounded by Cg ^ nia.x{\x^ys\ + a\y^\'^) , 
with the maximum taken over the same region. The free parameters a and b are 
used to control the Lipschitz constants Lf and Lg. Here A=l, implying ^^=£. 
Therefore holds for 0<5t<e and dm holds for e<dt< 2e. 

We first investigate how \E2\ scales with the macrostep At. We ensure that the term 
proportional to C*t^ is the dominant term in (|5.ip by initializing the fast variables 
on the approximate slow manifold with |d'^''^|=0. Figure ISTT] illustrates the linear 
scaling of \E2\ with At. System parameters are a — 1, 6 = 0.1. The scale separation 
parameter is e = 10^^. We used M — 90 microsteps with microstep size 6t — 0.1e and 
St =1.6£. The number of iterations n varied from 48 to 918 to keep T = nt^ = l 
fixed for all values of At. Initial conditions are chosen to lie on the approximate 
slow centre manifold with y° = l, = sin^(O.l). The Lipschitz constants are Lg = 2 
and Lf = 0.2, the bound on the vector field of the slow dynamics is Cg = l, and the 
maximal second derivative of the reduced slow dynamics is C* = 2. 

We now present results for the scaling of \E2\ with the time scale parameter e. To 
focus on the linear scaling suggested by the term LgLfCgS in (j5.1|) . we will have to 
control the term proportional to and the C*tA term. The |(i" |-term is controlled 
by setting initial conditions on the centre manifold and employing M^l, allowing 
for relaxation to the centre manifold. The C*t/\ term we control by choosing At<e, 
violating assumption (A6) (note that this implies 5t<£). System parameters are 
a = 0.1, b=l. Initial conditions are y*' = 5, a;" = sin^(5), i.e. |(i°'*'|=0. We used 
A/ = 100 microsteps with microstep size 5t=\0^^, macrostep size At ^10^^, for 
71 = 50 iterations of PI with T = 0.01. The Lipschitz constant for the slow dynamics is 
Lg = 2 and for the centre manifold is Lf = l, the bound on the vector field of the slow 
dynamics is Cg = 7, and the maximal second derivative of the reduced slow dynam- 
ics is C* = 6. Figure [12] illustrates the linear dependence of \E2 \ on e in this situation. 
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-5 -4 -3 

log(At) 



Fig. 5.1: Plot of logl^Jl versus logAt for fixed time of integration T=l. The dots 
represent results from numerical PI simulations of the system (|5.3|) - ()5.4p . with the 
crosses representing results from a system with St = 0.1e and the circles representing 
a system with St = 1.6e. The dashed lines are linear regression lines with a slopes of 
1.02 and 1.07, respectively. 



-11 



X 



X 



-9 -8 -7 

log(e) 

Fig. 5.2: Plot of logji?^*! versus loge. The dots represent results from numerical PI 
simulations of the system ()5.3|) - (|5.4p . The line is a linear regression with a slope of 
0.947. 



We now illustrate the linear scaling of \E'2\ with the maximal distance of 
the fast variable from the approximate centre manifold after a macrostcp. To 
ensure that the error is not dominated by the initial initialization error Id"'*^], 
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we choose parameters which render the scheme unstable and violate assumption 
{A8)^ allowing for divergence of the fast variables from the centre manifold over 
the macrosteps, i.e. > Figure shows clearly the linear dependence 

of I E"^ I on I d" I . System parameters are a = 1 , 6=1. Initial conditions are = 1 , 
a;°e [sin2(l) + 0.01,sin^(l) + 0.5]. The scale separation is e = 10"^. We used M = 100 
microsteps with microstep size St = 0.01e and St = 1.99e, and n = 5 macrosteps 
with macrostep size At=10~^ implying T = 0.0055 or T = 0.1. The Lipschitz 
constants are Lg = l and Lf^l, the bound on the vector field of the slow dynamics 
is Cg = 290, and the maximal second derivative of the reduced slow dynamics is C* = 6. 



O 



-5- 



-X','' 



O 



.x 



_"7l \ \ \ \ \ I 

1.5 2 2.5 3 3.5 4 

log I I 

Fig. 5.3: Plot of log|i?^'| versus log|d"|. The dots represent results from numerical PI 
simulations of the system (|5.3p - (|5.4p . The crosses are results from simulations with 
(5i = 0.0l£ and the circles are results from simulations with (5i = 1.99£. The dashed 
lines are linear regression lines with a slope of 1.00 and 1.03, respectively. 



6. Discussion We have established bounds on the error of a numerical approx- 
imation of the solution of a multiscale system by Projective Integration. The error 
contains terms stemming from the inherent error made by reducing the full dynamics 
on the centre manifold as well as errors specific to the numerical discretization. In 
particular, the order of the numerical scheme features, as well as errors due to inaccu- 
rately approximating the dynamics on the centre manifold. Although the constants 
involved in our error estimates are not optimal, the numerical simulations suggest 
that the scaling obtained is correct. 

In future work it is planned to use the analytical results obtained here as well as 
a trivial extension of our results to seamless HMM and the results obtained for the 
non-seamless version of HMM (see 0) to shed light on the important question in 
what circumstances one method or the other may lead to better performance. These 
methods exhibit different error bounds due to their different weights as well as due to 
differing reinitialization procedures for the fast variables. 
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